Lattice and surface effects in the out-of-equilibrium dynamics of the Hubbard model 
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We study, by means of the time-dependent Gutzwiller approximation, the out of equilibrium dy- 
namics of a half-filled Hubbard-Holstein model of correlated electrons interacting with local phonons. 
Inspired by pump-probe experiments, where intense light pulses selectively induce optical excitations 
that trigger a transient out-of-equilibrium dynamics, here we inject energy in the Hubbard bands 
by a non-equilibrium population of empty and doubly-occupied sites. We first consider the case of 
a global perturbation, acting over the whole sample, and find evidence of a mean-field dynamical 
transition where the lattice gets strongly distorted above a certain energy threshold, despite the 
weak strength of the electron-phonon coupling by comparison with the Hubbard repulsion. 
Next, we address a slab geometry for a correlated heterostructure and study the relaxation dynamics 
across the system when the perturbation acts locally on the first layer. While for weak deviations 
from equilibrium the excited surface is able to relax by transferring its excess energy to the bulk, 
for large deviations the excess energy stays instead concentrated into the surface layer. This self- 
trapping occurs both in the absence as well as in the presence of electron-phonon coupling. Phonons 
actually enforce the trapping by distorting at the surface. 

PACS numbers: 71.10.Fd, 71.30.+h, 78.47.-p 



I. INTRODUCTION 

The transient dynamical behavior of correlated mate- 
rials optically excited far from equilibrium is currently 
attracting growing interest, due to the impressive ad- 
vances in time-resolved spectroscopy with femtosecond 
resolution. By shining the sample with intense ultra-fast 
pulses (pump) , one can trigger non-equilibrium transient 
states, whose physical properties are then recorded by a 
second pulse arriving at fixed time delay (probe). The 
unique feature of these experimental techniques is to give 
access to dynamical information unavailable to conven- 
tional time-averaged frequency domain spectroscopies- 1 
In addition, when irradiation is sufficiently strong, one 
can even stabilize transient states with fundamentally 
different physical properties^ thus paving the way to a 
complete control of material properties by light 4 As cor- 
related electron systems are often on the verge of a Mott 
metal-to-insulator transition, this portends the utmost 
important possibility of optically manipulating their con- 
ducting properties on ultra-short time-scales^ 

Motivated by these achievements, the research activity 
on transient ultrafast dynamics in correlated electronic 
systems has rapidly grown in recent years£~— From a the- 
oretical perspective, one can expect a non trivial and rich 
transient dynamical behavior to emerge, reflecting the 
complex interplay between electrons, phonons, spins and 
orbital degrees of freedom that characterizes the phase di- 
agram of these materials. On a more fundamental level, 
the crucial question concerns whether these experiments 
could allow to explore novel metastable phases of corre- 
lated quantum matter that can only be accessed along 



non-thermal pathways. 

A common wisdom is that the effect of perturb- 
ing the system by an ultra-short laser pulse can be 
qualitatively accounted for by an effective-temperature 
description*^— Within this picture, the injected energy 
would turn first, on few femtoseconds, into heat for the 
electron sub-system only. At later times, picoseconds, 
the electronic heat is gradually transferred to the lattice 
so that, eventually, the whole system flows to a thermal 
state at higher temperature than the initial one. Un- 
der such a thermodynamic assumption, optical pumping 
should mimic the role of heating, hence allow accessing, 
possibly much faster, all phases that are reached upon 
raising temperature at equilibrium. In the specific case of 
correlated materials, this entails the possibility of photo- 
inducing metal-to-insulator transitions. Indeed, there ex- 
ist many examples of Mott insulators that can be driven 
metallic upon increasing temperature, like e.g. V20 o 14lls 
and VO^, and, vice versa, metals that turn Mott insu- 
lating upon heating, like the same V2O 3 14 ' 15 at higher 
temperatures, or like doped manganitesJ^ 

However, a deeper thought of what is known about 
correlated systems in equilibrium already raises questions 
on this point of view. Indeed, according to this picture 
one must conclude that energy-pumping, assumed to be 
equivalent to temperature raising, should make a metal 
less metallic and a band insulator less insulating. It is 
believed-^ that a correlated metal near a Mott transition 
actually shares properties of both metals and insulators; 
itinerant quasiparticles narrowly centered around the 
chemical potential coexisting with incoherent atomic-like 
high-energy excitations, the so-called Hubbard bands. If 
intense light exposure is the same as heating, and since 
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the light is selective via its frequency and polarization, 
then one could envisage the quasiparticles or else the 
Hubbard bands being heated first, which would corre- 
spond, respectively, to conductivity decrease or increase. 
Such a non-monotonous behavior would however contrast 
the effect of raising temperature at equilibrium, which 
is supposed to always lower conductivity^ The above 
observation thus challenges the picture of pump-probe 
experiments as effective thermodynamic perturbations. 

A different perspective, pointing toward an intrinsic 
kinetic nature of these experimental settings, is offered 
by the intense research activity around the non equi- 
librium dynamics of closed isolated many-body systems, 
which has recently attracted lot of interest in the differ- 
ent context of cold atoms trapped in optical lattices.— 
In this respect, it is by now well established that, when 
driven out of equilibrium by intense sudden perturba- 
tions, strongly correlated systems can be trapped into 
long-lived metastable states that differ qualitatively from 
their equilibrium counterpart. Example along this line 
is provided by the single band Hubbard model, likely 
the simplest model to describe strong correlation physics. 
Different theoretical approaches^— have shown, for ex- 
ample, that a sudden increase of the Hubbard repul- 
sion drives the system into a long-lived metastable state 
which, although highly excited, shows intrinsic features 
of a zero temperature metallic state, rather than in- 
coherent finite temperature effects as one would have 
guessed by thermodynamic arguments. Seemingly, sud- 
denly switching on a large Hubbard repulsion stabilizes 
metastable phases rich of energetically unfavorable dou- 
bly occupied sites, which are kinetically blocked 2 ^ and 
unable to decays or even to coherently propagated A 
qualitative picture of the crossovers or genuine dynami- 
cal transitions between different metastable states in the 
Hubbard model driven by sudden quantum quenches has 
been recently obtained using a time-dependent extension 
of the Gutzwiller approximation (t-GA) . 2 ' 25 While miss- 
ing important quantum fluctuations, which are crucial 
for the long-time dynamics, this approximate scheme has 
been shown to capture important qualitative features of 
the intermediate time evolution, which is actually of in- 
terest in the description of the pump-probe dynamics. 

In this work, we aim to elaborate further on this out 
of equilibrium perspective by including additional in- 
gredients that might play an important role in model- 
ing pump-probe experiments on actual correlated mate- 
rials. Firstly, we add phonons to the half-filled single 
band Hubbard model and study the transient dynamics 
induced by a sudden perturbation. It is worth notic- 
ing that lattice vibrations play a crucial role in actual 
experiments by triggering selective perturbation for the 
electronic subsystem 26 , and their role in ultrafast pump 
probe experiments is a subject of current experimental 
interest.— ~— Here, we consider Einstein phonons coupled 
to the local charge and study the dynamics of the result- 
ing Hubbard- Holstein model using a suitable extension of 
t-GA. Although extremely simplified, this model repre- 



sents a first attempt to figure out how highly excited elec- 
trons succeed in transferring their excess energy to the 
lattice. Results reveal, akin to the pure Hubbard model, 
the existence of a metastable state for high enough exci- 
tation, where phonons get strongly displaced in spite the 
large Coulomb repulsion and in striking contrast to what 
one would have guessed in equilibrium. 

A second important ingredient that we add to the de- 
scription builds on the observation, recently reported 
in a number of theoretical investigations, that non- 
thermal metastable states are extremely sensitive to spa- 
tial fluctuations and prone to spontaneous generation 
of inhomogeneitieS) 23 ' 30 which could play a crucial role 
in the dynamics, in particular around dynamical transi- 
tion pointSi 25 i 31 To investigate this issue, we consider the 
same Hubbard-Holstein model at half-filling but now in a 
slab geometry that lacks translational symmetry, model- 
ing ultra-fast dynamics in correlated heterostructuresj 2 ^ 
We assume that, initially, only the surface layer is driven 
out of equilibrium and study by t-GA how the excess en- 
ergy is redistributed inside the bulk. Remarkably, if the 
energy initially stored on the surface exceeds a critical 
threshold, it remains trapped on the uppermost layers. 
Concomitantly, the phonons distort at the surface, pro- 
viding a further trapping potential. This result demon- 
strates not only the importance of inhomogeneities, but 
also suggests that, under specific circumstances, the lat- 
tice might not provide a dissipative bath to speed up 
relaxation, but rather play the opposite game to slow 
down thermalization. 

The paper is organized as follows. In section [TT] we in- 
troduce the model and an out-of-equilibrium version of 
the Gutzwiller approximation that may cope with the 
electron-phonon coupling. In section III Al we show that 
the method is equivalent to the mean-field approximation 
applied to a model of free electrons coupled to phonons 
and Ising spins. In section Mil we move to discuss the 
results for two different cases. First, in section UlI Al we 
study the time-evolution when the whole bulk is suddenly 
driven out-of-equilibrium. Next, in section Till Bl we con- 
sider the situation in which an external pulse only excites 
the surface layer, and study if and how the surface can 
relax by transferring energy to the bulk. Finally, section 
IIVI is devoted to concluding remarks. 



II. THE MODEL AND THE GUTZWILLER 
APPROXIMATION 

We consider a half-filled Hubbard-Holstein model de- 
scribed by the Hamiltonian 

h = -t E kW+^O + ^E^- 1 ) 2 

<ij><j i 

+|E {p1 + x D -»E Xi fa- 1 )' 

i i 
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where c\ a (c ia ) creates(annihilates) an electron with spin 
a at site i, x$ is the phonon displacement at that site 
and pi its conjugate variable. The hopping is restricted 
to nearest neighbors and rii is the electron number op- 
erator. We note that ([T]) is invariant under particle-hole 
transformation provided Xi — > —Xi. 

In the following, we study the unitary dynamics in- 
duced by the Hamiltonian (p} using the Gutzwiller vari- 
ational scheme introduced at equilibrium by Barone et 
al— and extended to the time dependent case following 
Ref. H3- We emphasize at this point that, in real experi- 
ments, the light pulse couples via the vector potential to 
the electronic degrees of freedom. While in principle the 
variational description could be extended to include this 
feature, here we assume for the sake of simplicity that 
the effect of the pump is mainly to induce an initial non 
equilibrium distribution of electronic degrees of freedom, 
whose dynamics is then driven by the Hubbard-Holstcin 
Hamiltonian. Furthermore, while in real experimental 
settings the system is always in contact with a thermostat 
that eventually allows the injected energy to flow away, 
here we assume the whole system made by electrons and 
lattice to be isolated. While in different contexts, e.g. 
when current-carrying stationary states driven by static 
electric fields are present, this assumption may be highly 
questionable, here we stress that our focus concerns the 
transient relaxation dynamics on time scales of electronic 
and phononic degrees of freedom. As the coupling with 
the environment is typically very weak, we do not be- 
lieve that this assumption can qualitatively change the 
physical picture that we will draw. 

With these assumptions on the theoretical side, we 
now introduce our time-dependent variational wave func- 
tion for the dynamics of the Hubbard- Holstein model. ([1]) 
Specifically we write 

i*(i))=n>(^*) i*o(*)>, a) 

i 

where | ^o(t)) is a time-dependent Slater determinant, 
to be determined variationally, and Vi(xi,t) a time- 
dependent electron operator at site i that depends explic- 
itly on the phonon coordinate Xi. We define, neglecting 
the index i for convenience, 

v = V2Mx,t) |o)<o|+V2 0x(M)(rt)<t| + |+>Ul) 

+V2<f> 2 (x,t) |2)(2|, (2) 

where 4> n {x, t) are site-dependent phonon wave- functions, 
and | r) (r | is the projector onto the site being empty, 
r = 0, singly-occupied by a spin up, T =^, or down, 
r =1, electron, or, finally, doubly-occupied, T = 2. 
Particle-hole symmetry implies that, under n — >■ 2 — n, 
4> n (x,t) 4>2- n (-x,t), namely 

<po(x,t) = <j} 2 {-x,t), 
<j>i(x,t) = (j)i{-X,t). 



We evaluate average values on the wave-function (fTJ) 
by means of the Gutzwiller approximation, following 
Refs. [H and HH, which amounts to impose that 

J dx \^{ Xl t)\ 2 + \ct) 1 {x,t)\ 2 = 1. 

The above condition implies that the average over the 
Slater determinant | ^o(t)} and the phonons of the oper- 
ator that remains after extracting from Vi(xi,tyVi(xi, t) 
any two fermionic operators vanishes identically. This 
property allows to evaluate explicitly all average values 
on the wavefunction (fTJ) in the limit of infinite lattice- 
coordination ) 33 ' 34 although it is common to use the same 
results also for lattices with finite coordination numbers, 
hence the name Gutzwiller approximation. 

Within the Gutzwiller approximation, the average 
value of the Hamiltonian (p} on the wavefunction | Vt'(i)) 
can be shown to coincide with the average on | ^o(f)) of 
the Hamiltonian 

<ij>(T 

+ dx(u + 2gx) \<j> oi (x,t)f (3) 

+ -^^2 / dx4> ni (x,t)* h{x) 4> m {x,t) , 

i n=0,l 

where h(x) = ( — d^. + x 2 ). The parameters 

Ri{t) = J dx(<pu{x,t)*(poi(x,t) +c.c), (4) 

are commonly interpreted as the amplitudes of quasi- 
particles at sites i, hence if* as their effective non- 
interacting Hamiltonian with renormalized hopping 

The variational principle that we assume is the saddle 
point of the action S = J dt£(t), i.e. SS — 0, whose 
Lagrangian is 

£(t)=«<*(t) | *(«)>-<*(*) | H | *(<)), (5) 

which, within the Gutzwiller approximation^ reads sim- 
ply 

£(t) = X! / dx(f> ni (x,t)* 4> ni {x,t) (6) 

i n=0,l 

+i(Mt) I *o(*)> - (*o(t) I H,(t) | *„(*)). 

We define on each site i a normalized two-component 
spinor 

so that 



4 



and further introduce Pauli matrices er a , a = x,y,z 7 
which act on the two components of the spinor. With 
the above notations, the saddle point equations read: 

n.n. of i 



i I *o> 
where 



h (x) | $<) -t J2 R i W *J ° X I $ * 



Ri 



+ -(U + 2gx) (l-a z ) | $*), 
H* | *o), 



(*< I <? I *i>, 

E (*o i + H - c -) i *o) 



(8) 
(9) 

(10) 
(11) 



It is worth noticing here that, when the electron-phonon 
interaction vanishes, the two subsystems decouple and 
the above dynamics reduces, for the electronic degrees 
of freedom, to the one studied in Ref. [HI for the simple 
Hubbard model. In the general case, one has to inte- 
grate the equations of motion starting from initial values 
for the spinor wave- functions and the Slater determinant. 
In order to integrate the spinor part, we follow the ap- 
proach outlined in Ref. [3^ and project each component on 
the basis of eigenfunctions of the harmonic oscillator, the 
Hcrmitc functions <p n (x), namely we write for v = 0, 1 



MM) = E C n{ t ) l Pn{.x) 



(12) 



and obtain time dependent equations for the complex co- 
efficients c v n (t) by plugging this expansion into equation 
©. In practice, we truncate the basis set to a finite 
number of coefficients n = 0, . . . , iV& and check that con- 
vergence is guaranteed by choosing N b ~ 60. All cal- 
culations that are presented here have been performed 
on a cubic lattice using U = 12i and the phonon fre- 
quency u) — t. An important scale of energy is the value 
of the critical U c at the equilibrium Mott transition in 
the absence of phonons. In the cubic lattice and within 
the Gutzwiller approximation U c = 16i, whose inverse 
we shall use as the unit of time. As we mentioned, the 
Eqs. (|5j)- (fTTj) are strictly valid only in lattices with infi- 
nite coordination numbers, therefore our use in a cubic 
lattice is just an approximation. 



A. The Gutzwiller approximation as a mean-field 
theory 

The Eqs. © and §E§ resemble time-dependent mean- 
field equations, with the Schrcedinger-like evolution of 
| \l/o) that depends implicitly on the average values of 
selected operators over the wavefunctions | <&,), and vice 
versa for the latter ones. Indeed, one recognizes readily 
that, given the Hamiltonian 



Hi 



E 

<ij>cr 



(4. 



iE0+ 2 n) (w) 



(13) 



where erf, a = x,y,x, are Ising variables defined on each 
site i, and assuming a factorized wave-function 



where 



=| Ising+phonons) x | electrons) 



Ising+phonons) = | Ising+phonons) 



(14) 



the same variational principle 8S — that we applied 
before would lead right to Eqs. (O and ©. The Hamil- 
tonian (TlBf thus extends to the Hubbard- Holstein model 
the mapping derived in Ref. [H for the simple Hubbard 
model. We just recall that the mapping states that, if 
Z is the partition function of the original model with 
the Hamiltonian H of Eq. (p}, and Zj that one of the 
Hamiltonian Hi of Eq. (|13p , then, in the limit of infinite 
coordination lattices and at particle-hole symmetry, 



JV 



(15) 



where N is the number of sites^ Essentially, the map- 
ping demonstrates that the constraint required to imple- 
ment the so-called slave-spin representation of the Hub- 
bard modefc^~— is actually unessential in the limit of in- 
finite lattice-coordination and at particle-hole symmetry. 
The advantage of dealing with Hj instead of the original 
Hamiltonian is that it provides a simple framework to 
disentangle already at the mean field level the quasipar- 
ticle degrees of freedom, the fermionic operators, from 
the Hubbard bands, the Ising variables. 

The chosen factorization (|14[) . where the phonon de- 
grees of freedom are entangled with the Hubbard bands 
and both influence in a mean-field fashion the quasi- 
particles, is actually inspired by the DMFT result that, 
for large repulsion and weak electron-phonon coupling, 
phonon signatures are hardly visible in the quasiparti- 
cle spectrum but quite evident in the Hubbard bands.^- 
Different choices could be more appropriate in different 
contexts or easier to deal with, as the extreme factoriza- 
tion | =| Ising) x | phonons) x | electrons). 



III. RESULTS 

We shall now analyze the time-dependent mean field 
equations and (|S| that describe within t-GA the evo- 
lution of a variational wave-function under the action of 
the Hubbard-Holstein Hamiltonian, or equivalently the 
Hamiltonian (fl"3"|) . 

We will assume that the pump that drives the system 
out-of-equilibrium is selective in the sense that it only 
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injects energy in the Hubbard bands, i.e. in the Ising 
subsystem, specifically increasing the concentration of 
doubly-occupied sites (doublons), hence of empty sites 
(holons) because of particle conservation. In the Ising 
language, it corresponds to assuming that initially the 
average values of of are lower than those at equilibrium. 
We note that the equilibrium conditions are obtained by 
replacing the time-dependent mean-field equations Q 
and {§1 with stationary mean-field equations. In partic- 
ular, the equilibrium values of | <J>i) are the lowest energy 
eigenstates of the right hand side of Eq. © , which must 
be self-consistcntly determined since the effective Hamil- 
tonian depends on Rj — ($, | a x \ 

As we mentioned earlier, in real experiments the light 
pulse couples via the vector potential to both Hubbard 
bands and quasiparticles. Therefore the above assump- 
tion is only an approximation, whose validity we intend 
to weight up in the near future, while, in the present 
work, we shall keep assuming that the initial state is just 
characterized by an out-of-equilibrium equal population 
of doublons and holons. We will consider first the case 
in which such a population is uniformly distributed over 
the whole sample, and next move to inhomogeneous sit- 
uations. 



A. Whole bulk driven out-of-equilibrium 

Let us therefore consider the Hubbard-Holstein Hamil- 
tonian (TTJ) at half-filling and assume that the system is 
initially prepared with a uniform concentration of dou- 
blons and holons higher than at equilibrium. The system 
is then let evolve, its time evolution being approximated 
within t-GA by Eqs. (|9]) and (©. This case is actually 
similar to the quench described in Ref. l24t the new in- 
gredient being just the electron-phonon coupling. 

We first need to determine the equilibrium condition in 
the presence of the electron-phonon coupling. This is ac- 
complished by a self-consistent iterative mapping method 
similar to the one described in Ref . H(J The outcome is a 
homogeneous wave-function, with the same spinor at any 
site, and the Slater determinant that is just the uniform 
ground state of the hopping energy. The concentrations 
of doublons and holons are then artificially augmented 
by the same amount while keeping the phonon wave- 
function unaltered. This is accomplished by the scale 
transformation <fio(x) — » Xo4>o(x) and 4>i(x) — > Xi4>i(x), 
with Ao > 1 and Ai < 1 such that normalization is main- 
tained, 

J dx\l\Mx) | 2 +A 2 | </>iO) | a =l, 

but the concentration of doubly occupied and empty sites 
is increased. The system is then allowed to evolve as pre- 
viously explained. The novelty with respect to Ref. 
is that we can now monitor how the energy, initially in- 
jected in the electron subsystem only, is transferred to 
phonons. We note that, because translational symmetry 




t (u-») 

FIG. 1: The time evolution of the parameter Ri = R, Eq. (|l(jp 
for two different concentrations SD of injected doublons, one 
below and the other above the critical point, see Fig. [2] Below 
the critical point, SD — 4% (blue curve), R oscillates around a 
finite value, while, above, SD = 17% (red curve), it oscillates 
between +1 and -1 with zero average. 

is preserved by the time evolution, the effective Hamil- 
tonian H*(t) in Eq. © has R^t) = Rj(t) = R(t), V 
hence describes at all times a simple tight-binding model 
with uniform time-dependent nearest neighbor hopping. 
As a result, the Slater determinant that is initially the 
lowest energy eigenstate of the hopping, does not change 
in time, hence cannot provide dissipative channels for the 
spinor evolution. For this reason, the dynamics of both 
electronic and phononic observables lack relaxation to a 
steady state but rather shows undamped coherent oscil- 
lations. Still, as shown in Ref. the mean field dynam- 
ics captures important features of the non-equilibrium 
problem and provides a qualitatively correct picture of 
the short-to-intermediate time dynamics. 

In Fig. [T] we plot the time evolution of the renormal- 
ization factor Ri(t) = R(t), Eq. (TTU)) . which shows two 
distinct regimes of oscillations depending on the amount 
of doublons injected, SD, which measures the strength of 
the non equilibrium perturbation and that we define as 
SD = D(t = 0) - L>cq., with D(t = 0) the initial value 
and -D e q. the equilibrium one. 

For small perturbations, R oscillates around a finite 
average, while, upon increasing SD above a threshold, it 
oscillates from -1 to +1, with average zero. In the Ising 
model language of section Hi A[ this behavior is represen- 
tative of the transition from the ordered phase, (a x ) ^ 0, 
to the disordered one, (a x ) = 0. A similar dynamical 
transition was observed in Ref. [U for the pure Hubbard 
model without electron-phonon coupling. The mean field 
coherent oscillations, although artificial, reflect the real 
tendency of the system to be trapped into long lived pre- 
thermal metastable states, whose properties are correctly 
captured by the long-time averages of the mean field dy- 
namics. 

With this insight, we plot in Fig. [5] the time- averaged 
double occupancy as a function of the concentration of 
injected doublons SD. The first observation is that the 
electron-phonon interaction does not change qualitatively 
the behavior with respect to the Hubbard model alone, 
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FIG. 2: The time average value of the percentage of dou- 
bly occupied sites as function of the percentage SD of dou- 
blons, hence holons, injected in the initial state with respect 
to the equilibrium value. We plot three values of the electron- 
phonon couplings, g = (black), g = O.lt (blue) and g = 0.2t 
(red). Larger g correspond to more pronounced kinks near 
the critical point, which is identified by the point at which 
the double occupancy drops down. 

see Ref. HH namely, we still find two distinct regimes sep- 
arated by a critical point where the double occupancy 
goes to zero, although numerically we cannot hit the 
precise value when this occurs. We note that the lo- 
cation of the critical point is not appreciably affected 
by phonons because of the tiny electron-phonon coupling 
(g 2 /u ~ 10~ 3 U). 

We find that the transition occurs right when the ini- 
tial energy happens to coincide with the equilibrium en- 
ergy at the Mott transition ; 24 ! 31 which is simply the zero 
point energy of the phonons for our model Hamiltonian 
([T]) and within the Gutz wilier approximation. In fact, 
one readily realizes that Eqs. ((U) and © admit another 
stationary point besides the one that corresponds to the 
equilibrium condition, namely <fio{%) = hence R = 0, 
with energy just w/2. We finally mention that, unlike in 
the absence of phonons^ here the critical point is not 
associated to an exponential relaxation towards a station- 
ary state that seems to be a characteristic of integrable 
dynamics^ which is presumably not our case. 

We now move our attention to the phonon sector, in 
order to unveil the entanglement between the electrons 
and the lattice as the former are driven out of equilib- 
rium. A natural quantity to look at would be the average 
lattice displacement (xi), which however is constrained to 
be zero on average by particle-hole symmetry. Still, we 
can define as a measure of the effective displacement the 
average of the operator q{ = Xi(rii — 1), which is just the 
electron-phonon coupling operator. In Fig. [3] we plot the 
relative variation of the time-averaged effective displace- 
ment 

q* = lim - / dt <?<(*)>, (16) 
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FIG. 3: The time average value of the lattice distortion, de- 
fined as the relative variation with respect to the equilibrium 
value, as function of the percentage of doublons 8D injected 
in the initial state with respect to the equilibrium value. The 
curve that is more singular near the critical point corresponds 
to g = O.lt (blue), the other to g = 0.2t (red). 

i.e. (g* — 9cq.)/<7cq., where q cq , is the equilibrium value, as 
function of the concentration of injected doublons. We 
note that, at small concentrations, the displacement is 
mostly unchanged from its equilibrium value. However, 
for higher concentrations past the critical point, the dis- 
placement starts increasing substantially; a growing dis- 
tortion being a way to store the initial excess energy. 




injected doublons (%) 



FIG. 4: Color plot of the Fourier transform (in arbitrary 
units) of the time evolutions of the double occupancy and the 
effective phonon distortion, as function of the concentration 
of injected doublons and of the frequency in units of U c . 

Although the gross behavior seems not to be affected 
by phonons, there are details which feature their pres- 
ence. In particular, we note some anomalies, tinier in 
Fig. [2] and more visible in |3J These anomalies appear 
when the oscillation frequency of the electronic dynam- 
ics, which decreases on approaching the critical point, 



7 




FIG. 5: Slab geometry for simulating the hypothetical exper- 
iment in which only the surface layer is driven out of equilib- 
rium 



hits the renormalized phonon frequency, or a multiple 
of it. Since the latter is small, these resonances occur 
near the critical point. This is evident in Fig 01 where 
we draw by a color plot the spectral decomposition of 
the time evolutions both of the double occupancy and 
the phonon distortion as function of the concentration of 
injected doublons and of the frequency of the signal. In 
particular, we observe the avoided crossing between the 
two lowest frequencies around SD ~ 5% that causes the 
kink visible in FigJ3] Among these two frequencies, the 
lowest one is visible mostly in the dynamics of q(t), hence 
can be regarded as a renormalized phonon frequency that 
blue-shifts upon increasing SD 7 i.e. the energy injected 
into the system. 



B. Surface driven out-of-equilibrium 

Let us now consider a slab geometry as depicted in 
Fig. [SJ denoting by z the direction perpendicular to the 
surface, which lies in the xy-p\ane. This setting allows 
us to mimic the non-equilibrium dynamics across a cor- 
related heterostructure, that has recently attracted ex- 
perimental interest^ We consider a system described by 
the Hubbard- Holstein model and consider a perturbation 
acting only at the surface layer, by triggering an out- 
of-equilibrium population of doublons and holons, while 
keeping the bulk in its equilibrium ground-state configu- 
ration. This initial state is then let evolve and its time- 
evolution is approximated by the Eqs. ([9]) and (JU). 

This particular geometry has the additional complica- 
tion that, at equilibrium, the optimized | $j) are layer 
dependent and the optimized Slater determinant is not 
uniform anymore. Therefore, the first step we need to 
undertake is solving the equilibrium problem, which we 
accomplish by the method developed in Ref. HO- Because 
of the slab geometry, we can choose a basis of single- 
particle wave-functions for building the Slater determi- 
nant defined by 

e zk-r 

V>ek(r,i) = —j= 1pek(t), 
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i i i i i i i i i i i i i i i i i i i i i 

50 100 150 200 

time (U; 1 ) 

FIG. 6: Time evolution of the percentage of doubly occupied 
sites on three different layers, indicated in the figure, for small 
percentage of injected doublons. The time is measured in 
units of the inverse of U c . The simulation is performed with 
a slab of 100 layers. 

where r is the space coordinate and k = (k x ,k y ) the 
momentum in the xy-plane, which is assumed to contain 
A lattice sites, while i = 1,...,N is the layer index, 
and typically we used N = 100 layers. Since there is 
translational symmetry in the xy-plane, we can choose 
the spinor | <F;) to depend only on the layer index i. 
Then, the stationary solution of ([8]) and (O amounts to 
solve at fixed Ri the eigenvalue problem 

e4>ek(i) = tRiei i ifj e u(i)~tR i ^ Ri+aipck{i + a)-, 

a=±l 

(17) 

where e k = — 2 (cos 

kx ~t" cosfcy), with the boundary con- 
dition Vek(O) = ipek(N + 1) = 0. The lowest energy 
eigenfunctions e < cp, cf — because of particle- hole 
symmetry, are then used to define the Slater determi- 
nant | ^o) and the average hopping between layer i and 
i + a 

Wj-^i+a = \ ^ ^2 (^ek(i)*ipek(i + a) + c.c) , (18) 

e<0 k3c<0 

as well as the average hopping within layer i 

W ^ = jJL E * KekWI 2 . (19) 

e<0 k9e<0 

These parameters are used to solve the spinor eigenvalue 
problem 

E | = | h{x) | $i) - tRi w^i a x | *<) 

-t ^ R i+a Wi^i+a ^ | ^i+a) 
a=±l 

+ l(U + 2gx)(l~a z ) \ <D 4 ), (20) 
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FIG. 7: Same as Fig. [5] at higher percentage of injected dou- 
blons. 
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FIG. 8: The lattice distortion on three different layers with 
the same amount of injected doublons as in Fig. [7J 



whose lowest energy solution defines new parameters 
Ri = ($j | a x | $j) that are used to solve again Eq. ([17]) 
and so on, until convergence is reached ! 39 ' 

The breaking of translational symmetry by the pres- 
ence of the surfaces is actually amplified by electron cor- 
relations that create a surface dead iayer ) 40 i 41 with sup- 
pressed double occupancy hence reduced hopping renor- 
malization parameters Ri. The dead layer penetrates 
inside the bulk over a length proportional to the Mott- 
transition correlation length i 40 ' 41 

Given this starting state, we suddenly increase the pop- 
ulation of doublons and holons on the first layer i = 1 
and let the system evolve. Essentially, we simply turn 
Eqs. (fT7|) and ([20]) into self-consistent non-linear time- 
dependent Schrcedinger equations that we solve numer- 
ically. Unlike in the homogeneous case of section 1111 Al 
here the Slater determinant evolves with time because 
the wave- functions ip e k(i, t) acquire a non-trivial time de- 
pendence, which provides additional dissipative channels 
that were previously absent. In other words, the hopping 
parameters Wi-yi+ a (t) and Wi^i{t) defined in Eqs. ((18"]) 
and (I19[) become time-dependent and influence the evo- 
lution of | see Eq. (|20p. which in turns affects il>ek{i, t) 
via the parameters Ri(t). 

The mutual feedback between | <F ) and the | 3>j)'s 
brings about a non-trivial dynamics much richer than 
in the example discussed in section IIII Al Nevertheless, 
even in this case we do find two completely different dy- 
namical behaviors, depending on the amount of injected 
doublons. For small values, the perturbed surface layer 
is able to relax by dissipating its excess energy inside 
the bulk, see Fig. [6] Indeed, the time-average values of 
the double occupancies on each layer tend towards their 
equilibrium values, which, as we mentioned, are lower the 
closer the layer to the surface. We emphasize that here, 
as opposite to the previous case and to the case of global 
quantum quenches, the perturbation is local and the en- 
ergy injected does not scale with the system size. As 
a result, the relaxation dynamics we find in this regime 



is toward the equilibrium ground state and no heating 
or finite temperature effects are expected in the long- 
time limit. This is a specific example of a local quantum 
quench and shows that our time dependent Gutzwiller 
approximation, with the above mentioned feedback be- 
tween variational parameters and Slater determinant, is 
able to describe thermalization. We also mention that, 
working with a finite-size geometry, recurrence effects 
are present at long enough times, when the perturba- 
tion reaches the opposite surface and starts oscillating 
back and forth. We expect that by taking the thermody- 
namic limit in the z direction these finite-size effects will 
be washed away. Nevertheless, even for finite lengths the 
relaxation and the trend towards equilibrium are clearly 
evident. 

Upon increasing the concentration of injected doublons 
a different dynamical behavior emerges. In particular, 
above a certain threshold the excitation remains trapped 
near the surface, see Fig. [7J This is quite remarkable 
because, as we mentioned, the Slater determinant now 
adjusts to the spinors | $i) during the time evolution, 
hence could in principle absorb the excess energy and 
transfer it in the interior of the bulk. What actually 
happens is that the parameters Ri and -Rj+i of adjacent 
layers interfere destructively, i.e. oscillate out-of-phase, 
at some i near the surface, effectively suppressing the 
quasiparticle inter-layer hopping t ii+ i = Ri Ri+x t, hence 
cutting layer i from the rest of the bulk. This anomalous 
trapping exists also in the absence of electron-phonon 
coupling, hence it is primarily an electronic effect, pre- 
sumably the dynamical counterpart of the surface dead 
iayer at equilibrium ^ 40 ' 41 What changes at finite electron- 
phonon coupling is that this phenomenon is accompanied 
by a lattice deformation, also localized on the uppermost 
layers, see Fig. [5] 

The physical picture that emerges can be visualized 
much better by the long-time layer-dependent profiles of 
the percentages of doubly occupied sites and of the dis- 
tortion, shown in Fig. [S] We observe that the deviations 
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FIG. 9: The percentage of doubly occupied sites and of the 
lattice distortion at each layer at very long times. 



of both quantities with respect to equilibrium are indeed 
concentrated just near the surface, while the bulk is prac- 
tically unaffected. Also instructive are the profiles of the 
intra-layer and inter-layer hopping renormalization fac- 
tors, Rf and RiRi + i, shown in Fig. [TUJ We note that 
the first layer has a much larger hopping renormalization 
factor, Rf, than at equilibrium, when it would be very 
small due to the dead layer phenomenon ! 40 ' 41 However, 
this layer is practically decoupled from the second layer, 
i?i i?2 being vanishingly small. 

We end by mentioning that, in contrast to the case of 
section IIII A[ the two different regimes that we observe 
seem not to be separated by a genuine dynamical criti- 
cal point, but rather by the dynamical counterpart of a 
first order phase transition. Indeed, in the intermediate 
regime the system does not show a well defined behav- 
ior but instead oscillates between the two distinct phases 
above. 



IV. CONCLUSIONS 

In this work we have studied the real time dynamics 
of the Hubbard-Holstein model at half-filling by a very 
simple extension of the Gutzwiller approximation in two 
different toy cases: (i) a bulk system is prepared with an 
equal out-of-equilibrium population of doubly occupied 
and empty sites and let evolve in time; (ii) a slab is con- 
sidered and it is assumed that only the surface layer is 
initially driven out-of-equilibrium. 

In case (i) we find similar results as in the quantum 
quench of the pure Hubbard model: a dynamical criti- 
cal point that separates two different regimes. The novel 
feature introduced by the phonons is the presence of a 
substantial phonon-displacement that occurs for large 
enough deviation from equilibrium, a remarkable out- 
come in that the electron-phonon coupling we consider 
is extremely small as compared with the Hubbard repul- 
sion. 

In the slab geometry (ii) we still find two dynamical 
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layer i 

FIG. 10: The intra-layer, bottom panel, and inter-layer, top 
panel, hopping renormalization factors at long times. 



behaviors, although this time without a true dynami- 
cal transition in between. If the energy injected at the 
surface is below a threshold, it is able to diffuse within 
the bulk and the system seems to relax towards the 
equilibrium ground state with the dead layer near the 
surface ] 40 ' 41 On the contrary, if the excess energy at the 
surface exceeds that threshold, it does not succeed any- 
more to diffuse in the bulk and remains concentrated 
practically at the surface, bringing about a substantial 
phonon displacement. Surprisingly, we finds that the 
first layer has a larger hopping renormalization factor 
than at equilibrium, which can be sustained because the 
layer effectively decouples from the rest of the system. 
However, we can not conclude that such an enhancement 
corresponds to an increased metallicity, which would be 
indeed a remarkable result. In fact, we tend to believe 
that the hopping renormalization as defined within the 
Gutzwiller approximation is a measure of the whole, co- 
herent plus incoherent, single-particle spectrum at low 
energy, not just of the quasiparticle coherent contribu- 
tion alone. Therefore, what we feel safe to state is just 
that low energy spectral weight grows in the first layer, 
whatever being its nature. 

We cannot exclude that such a long-lived localized ex- 
citation could indeed correspond to some kind of exciton 
already present in the equilibrium spectrum, which can 
be unveiled by our variational technique only because we 
are exploring the dynamics. It is also plausible that such 
a localized excitation exists just in correspondence with 
the surface dead iayer ; 40 ' 41 where the low-energy spectral 
weight is negligible hence there is room for excitons inside 
the preformed Mott-Hubbard gap. It is as well possible 
that our finding is actually related to the debated issue 
about the lifetime of doublons in the strongly interact- 
ing Hubbard model^ ' 22 ' 42 which could also be the clue to 
understand the lack of thermalization of highly excited 
states when correlation is strong. Further investigations 
with different and complementary approaches are needed 
to clarify these interesting issues. 
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